home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / fft.lha / fft / fft_output.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  5KB  |  183 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *                  Fast Fourier Transform
  6.  *
  7.  *        Return the transformation results in the form
  8.  *               the user wants them
  9.  *
  10.  ************************************************************************
  11.  */
  12.  
  13. #include "fft.h"
  14. #pragma implementation "fft.h"
  15.  
  16. #include <math.h>
  17. #include <Complex.h>
  18. #include <std.h>
  19.  
  20. /*
  21.  *-----------------------------------------------------------------------
  22.  *         Give a real/imaginaire part / absolute value
  23.  *            of the complex transform
  24.  */
  25.  
  26. void FFT::real(Vector& xf_re)        // [0:N-1] vector
  27. {
  28.   if( xf_re.q_lwb() != 0 || xf_re.q_upb() != N-1 )
  29.     _error("Vector [%d:%d] isn't fit for placing the resulting transform.\n",
  30.        "[%d:%d] vector was expected",
  31.        xf_re.q_lwb(), xf_re.q_upb(), 0, N-1);
  32.  
  33.   register int i;
  34.   register Complex * ap;
  35.   for(i=0,ap=A; i<N; i++)
  36.     xf_re(i) = (*ap++).real();
  37.   assert( ap == A_end );
  38. }
  39.  
  40. void FFT::imag(Vector& xf_im)        // [0:N-1] vector
  41. {
  42.   if( xf_im.q_lwb() != 0 || xf_im.q_upb() != N-1 )
  43.     _error("Vector [%d:%d] isn't fit for placing the resulting transform.\n",
  44.        "[%d:%d] vector was expected",
  45.        xf_im.q_lwb(), xf_im.q_upb(), 0, N-1);
  46.  
  47.   register int i;
  48.   register Complex * ap;
  49.   for(i=0,ap=A; i<N; i++)
  50.     xf_im(i) = (*ap++).imag();
  51.   assert( ap == A_end );
  52. }
  53.  
  54. void FFT::abs(Vector& xf_abs)        // [0:N-1] vector
  55. {
  56.   if( xf_abs.q_lwb() != 0 || xf_abs.q_upb() != N-1 )
  57.     _error("Vector [%d:%d] isn't fit for placing the resulting transform.\n",
  58.        "[%d:%d] vector was expected",
  59.        xf_abs.q_lwb(), xf_abs.q_upb(), 0, N-1);
  60.  
  61.   register int i;
  62.   register Complex * ap;
  63.   for(i=0,ap=A; i<N; i++)
  64.     xf_abs(i) = ::abs(*ap++);
  65.   assert( ap == A_end );
  66. }
  67.  
  68. /*
  69.  *-----------------------------------------------------------------------
  70.  *           Give only a half of the resulting transform
  71.  */
  72.  
  73.  
  74. void FFT::real_half(Vector& xf_re)    // [0:N/2-1] vector
  75. {
  76.   if( xf_re.q_lwb() != 0 || xf_re.q_upb() != N/2-1 )
  77.     _error("Vector [%d:%d] isn't fit for placing the resulting transform.\n",
  78.        "[%d:%d] vector was expected",
  79.        xf_re.q_lwb(), xf_re.q_upb(), 0, N/2-1);
  80.  
  81.   register int i;
  82.   register Complex * ap;
  83.   for(i=0,ap=A; i<N/2; i++)
  84.     xf_re(i) = (*ap++).real();
  85.   assert( ap + N/2 == A_end );
  86. }
  87.  
  88. void FFT::imag_half(Vector& xf_im)    // [0:N/2-1] vector
  89. {
  90.   if( xf_im.q_lwb() != 0 || xf_im.q_upb() != N/2-1 )
  91.     _error("Vector [%d:%d] isn't fit for placing the resulting transform.\n",
  92.        "[%d:%d] vector was expected",
  93.        xf_im.q_lwb(), xf_im.q_upb(), 0, N/2-1);
  94.  
  95.   register int i;
  96.   register Complex * ap;
  97.   for(i=0,ap=A; i<N/2; i++)
  98.     xf_im(i) = (*ap++).imag();
  99.   assert( ap + N/2 == A_end );
  100. }
  101.  
  102. void FFT::abs_half(Vector& xf_abs)    // [0:N/2-1] vector
  103. {
  104.   if( xf_abs.q_lwb() != 0 || xf_abs.q_upb() != N/2-1 )
  105.     _error("Vector [%d:%d] isn't fit for placing the resulting transform.\n",
  106.        "[%d:%d] vector was expected",
  107.        xf_abs.q_lwb(), xf_abs.q_upb(), 0, N/2-1);
  108.  
  109.   register int i;
  110.   register Complex * ap;
  111.   for(i=0,ap=A; i<N/2; i++)
  112.     xf_abs(i) = ::abs(*ap++);
  113.   assert( ap + N/2 == A_end );
  114. }
  115.  
  116. /*
  117.  *-----------------------------------------------------------------------
  118.  *        Perform sin/cos transforms of a real function
  119.  *              as a postprocessing of FFT
  120.  *
  121.  * Sine-transform:   F(k) = Integrate[ f(x) sin(kx) dx ], x = 0..Infinity
  122.  * Cosine-transform: F(k) = Integrate[ f(x) cos(kx) dx ], x = 0..Infinity
  123.  * Inverse
  124.  *   sin-transform:  f(x) = 2/pi Integrate[ F(k) sin(kx) dk ], k = 0..Infinity
  125.  *   cos-transform:  f(x) = 2/pi Integrate[ F(k) cos(kx) dk ], k = 0..Infinity
  126.  *
  127.  * Function f(x) is tabulated over the uniform grid xj = j*dr, j=0..n-1
  128.  * Function F(k) is tabulated over the uniform grid kj = j*dk, j=0..n-1
  129.  *                             n=N/2
  130.  * Source and destination arguments of the functions below may point to
  131.  * the same vector (in that case, transform is computed inplace)
  132.  */
  133.  
  134. void FFT::sin_transform(Vector& F, const Vector& f)
  135. {
  136.   are_compatible(F,f);
  137.   input_pad0(f);
  138.  
  139.   register int j;
  140.   register Complex * ap;
  141.   for(j=0,ap=A; j<N/2; j++)
  142.     F(j) = - (*ap++).imag() * dr;
  143.   assert( ap + N/2 == A_end );
  144. }
  145.  
  146. void FFT::cos_transform(Vector& F, const Vector& f)
  147. {
  148.   are_compatible(F,f);
  149.   input_pad0(f);
  150.  
  151.   register int j;
  152.   register Complex * ap;
  153.   for(j=0,ap=A; j<N/2; j++)
  154.     F(j) = (*ap++).real() * dr;
  155.   assert( ap + N/2 == A_end );
  156. }
  157.  
  158.  
  159. void FFT::sin_inv_transform(Vector& f, const Vector& F)
  160. {
  161.   are_compatible(F,f);
  162.   input_pad0(F);
  163.  
  164.   register int j;
  165.   register Complex * ap;
  166.   for(j=0,ap=A; j<N/2; j++)
  167.     f(j) = - (*ap++).imag() * 4/N/dr;    // 2/pi * dk = 2/pi * 2pi/N/dr
  168.   assert( ap + N/2 == A_end );
  169. }
  170.  
  171. void FFT::cos_inv_transform(Vector& f, const Vector& F)
  172. {
  173.   are_compatible(F,f);
  174.   input_pad0(F);
  175.  
  176.   register int j;
  177.   register Complex * ap;
  178.   for(j=0,ap=A; j<N/2; j++)
  179.     f(j) = (*ap++).real() * 4/N/dr;
  180.   assert( ap + N/2 == A_end );
  181. }
  182.  
  183.